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THE EFFICIENT EVALUATION OF THE HYPERGEOMETRIC 
FUNCTION OF A MATRIX ARGUMENT 

PLAMEN KOEV AND ALAN EDELMAN 



Abstract. We present new algorithms that efficiently approximate the hy- 
pergeometric function of a matrix argument through its expansion as a series 
of Jack functions. Our algorithms exploit the combinatorial properties of the 
Jack function, and have complexity that is only linear in the size of the matrix. 



1. Introduction 

The hypergeometric function of a matrix argument has a wide area of appli- 
cations in multivariate statistical analysis |17j . random matrix theory [7j, wireless 
communications [HI EL etc- Except in a few special cases, it can be expressed only 
as a series of multivariate homogeneous polynomials, called Jack functions. This 
series often converges very slowly ^3 El P- 390] , and the cost of the straightfor- 
ward evaluation of a single Jack function is exponential [3]. The hypergeometric 
function of a matrix argument has thus acquired a reputation of being notoriously 
difficult to approximate even in the simplest cases [21 EH ■ 

In this paper we present new algorithms for approximating the value of the 
hypergeometric function of a matrix argument. We exploit recursive combinatorial 
relationships between the Jack functions, which allow us to only update the value 
of a Jack function from other Jack functions computed earlier in the series. The 
savings in computational time are enormous; the resulting algorithm has complexity 
that is only linear in the size of the matrix argument. In the special case when the 
matrix argument is a multiple of the identity, the evaluation becomes even faster. 

We have made a MATLAB ^5] implementation of our algorithms available J3| • 
This implementation is very efficient (see performance results in Section H3), and 
has lead to new results |B] . 

The hypergeometric function of a matrix argument is defined as follows. Let 
p > and q > be integers, and let X be an n x n complex symmetric matrix with 
eigenvalues x\, X%, . . . , x n . Then 

(1.1) p it>K . . . , a p ; bl , . . . , b g ; X) = £ £ ^ V 1 >>> ' C ^ ^ 
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where a > is a parameter; n h k means k = (ki, K2, ■ . .) is a partition of fc (i.e., 
ki > K2 > . . . > are integers such that = Ki + K2 + ■ • • = fc); 

(1.2) (a)W= J] ( a -^ + ^'- 1 ) 

is the generalized Pochhammer symbol, and Cr (X) is the Jack function. 

The Jacfc function Ck(X) = (x\, X2, ■ ■ ■ ,x n ) is a symmetric, homogeneous 
polynomial of degree |k| in the eigenvalues x±,X2, ■ ■ ■ , x n of X Rem. 2, p. 228], 
|2()j . For example, when a = 1, (X) becomes the (normalized) Schur function, 
and for a = 2, the zonal polynomial. 1 There are several normalizations of the Jack 
function which are scalar multiples of one another: (X) is normalized so that 
X) K i-fc Cre^pO = (trX) fc ; in Scction|21we express (|1.1|) in terms of the Jack function 
Jk (X), which is normalized so that the coefficient of X1X2 • ■ -x\ K i is (|k|)!. The 
functions C^\X) and J^Hx) can be defined recursively, e.g., 

(1.3) JW(ii,u, • • • ,x n ) = Y, JfKxi, X2, ■ ■ ■ ,x n -i) ■ xW»\ ■ 0^, 

where /3 KM is a rational function of a (see Section [3] for details). The relationship 
(11.3(1 becomes key in achieving efficiency in our algorithms. 

The Jack functions C^\X) and Jk(X), and in turn the hypergeometric func- 
tion of a matrix argument, depend only on the eigenvalues x\,X2, ■ ■ ■ ,x n of X. 
Many authors, however, have found the matrix notation in Ql.lfl . and the use of a 
matrix argument to be more convenient. We follow the same practice. 

The hypergeometric function of a matrix argument is scalar-valued, which is 
a major distinction from other functions of a matrix argument (e.g., the matrix 
exponential), which are matrix-valued. The hypergeometric function of a matrix 
argument generalizes the classical hypergeometric function to which it reduces for 
n = 1. In general, however, there is no explicit relationship between these two 
functions for n > 2. 

We approximate the series l)l.l|) by computing its truncation for \k\ < m: 

m 1 ^( a ) / \( a ) 

(1.4) . . . , ap ; bl , . . . , b q ; X) = £ £ g^P^ • (X). 

fc=o K hfe 

The series (jl.lll converges for any X when p > q; it converges if maxj \xi\ < 1 
when p = q + 1, and diverges when p > q + 1, unless it terminates [171 p. 258]. 
When it converges, its K-term converges to zero as \k\ — * 00. In these cases (|1.4|) is 
a good approximation to 1(1. ljl for a large enough m. 

The computational difficulties in evaluating ((1.4(1 are: 

(A) the series (|1 . 1|) converges slowly in many cases ^B]; thus a rather large m 
may be needed before 1(1.4(1 becomes a good approximation to 1(1.1(1 : 

(B) the number of terms in 1(1.4(1 (i.e., the number of partitions \k\ < m) grows, 
roughly, as 0(e v ^™) (see Section |S}; 

(C) the straightforward evaluation of a single Jack function, Cl"'^) for \k\ — 
m, has complexity that grows as 0(n m ) [H]- 

^Some authors define the hypergeometric function of a matrix argument through the series 
11.11 for a = 1 |5] (4.1)] or a = 2 1171 p. 258] only. There is no reason for us to treat the different 
cr's separately (see also Q]|5]|7| for the uniform treatment of the different ct's in other settings). 
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While there is little we can do about (A) (which is also a major problem even in 
the univariate (n — 1) case |18|V or (B), our major contribution is in improving (C), 
the cost of evaluating the Jack function. We exploit the combinatorial properties of 
the Pochhammer symbol and the Jack function to only update the K-term in (|1.4() 
from the /z-terms, /i < k. As a result the complexity of our main algorithm for 
computing l|1.4|l . Algorithm 14. 21 is only linear in the size n of the matrix argument 
X, exponentially faster than the previous best algorithm (see Sections |21 El and 
16.31 for details). In the special case when X is a multiple of the identity, we present 
an even faster algorithm, Algorithm 14.11 whose complexity is independent of n. 

A number of interesting problems remain open. Among these are: 

• detecting convergence; 

• selecting the optimal value of m in l|1.4l) for a desired accuracy; 

• selecting the optimal truncation of the series Hl.lfl . 

We do not believe that a uniform answer to these problems exists for every a and 
every p and q. Therefore, we leave the choice of m and an appropriate truncation 
to the user. We elaborate more on these open problems in Section [3 

With minimal changes our algorithms can approximate the hypergeometric func- 
tion of two matrix arguments 



(1.5) vF^Hcu-vM.^X^Y) = > > 

v I p q y J-i" ' > Z-^i L^i ii/. uoj /. \\a) r". a >(T\ 

k=0 K hk K -{°1)k ■■■\°q)K G-K (1) 



(q) (l \(") r<(°>)i 



and more generally functions of the form 

oo 

(1.6) G(X) = ^^2a K C^(X), 



k=0 Khk 



for arbitrary coefficients a K at a similar computational cost (see, e.g., I|6.5f) in sub- 
section E2Jh 

In l|1.5fl and throughout this paper, we denote a vector (z 1; . . . , z t ) as z 1:t . 

This paper is organized as follows. We survey previous algorithms for computing 
the hypergeometric function of a matrix argument in Section [5] In Section [3] wc 
describe our approach in computing the truncation ljl.4|l . We present our new 
algorithms in Section and analyze their complexity in Section We present 
numerical experiments in Sectional Finally, we draw conclusions and present open 
problems in Section 



2. Previous Algorithms 

Butler and Wood used Laplace approximations to compute the integral rep- 
resentations [HI Thm. 7.4.2, p. 264]: 

^(a;c;X) = ^ f e^fdrtY)^ 

1 n n (C - a) Jo<Y<I 

x det{I -Yy-a-^idY), 
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valid for real symmetric X, SR(a) > 5R(c) > and 3?(c - a) > and 

2 if>M;c;X) = -/ det(I-XY)-» 

I }j (a)l n (c - a) Jo<Y</ 

x(detr) a - i * i det(/ - r) c - a -^(dr), 

valid for K(X) < /, 5R(a) > and 3?(c - a) > 

This approach, however, is restricted to the cases p = 1 or 2, q = 1, and a = 2. 

Gutierrez, Rodriguez, and Saez presented in ^U] (see also |T§] for the implemen- 
tation) an algorithm for computing the truncation (|1.4|) for a = 2 (then the Jack 
functions are called zonal polynomials). For every k = 1, 2, . . . , m, the authors form 
the upper triangular transition matrix |141 p. 99] K (indexed by all partitions of 
k) between the monomial symmetric functions (m K ) K \-k and the zonal polynomials 
(C K ) K \-k- Then for every partition k = (ki, K2, . . .) h k they compute 

(where /i ranges over all distinct permutations of k 2 , . . .) |22l p. 289]), and form 
the product (C K ) K |-/c = K ■ (TO K ) K i-fc- Computing the vector (m K ) K |_ m alone costs 
m { n+ m~ ) = 0{n m ) since every term in every m K is of degree m, and for every 
nonstrictly increasing sequence of m numbers from the set {1, 2, . . . , n} we obtain 
a distinct term in some m K . The overall cost is thus at least exponential (0(n m )), 
which explains the authors' observation: 

We spent about 8 days to obtain the 627 zonal polynomials of 
degree 20 with a 350 MHz Pentium II processor. 

In contrast, our Algorithm 14. 21 takes less than a hundredth of a second to do the 
same. Its complexity is only linear in n and subexponential in m (see Section EJ. 

3. Our Approach 

We make the evaluation of ™Fq (ai :p ; b\- q \ X) efficient by only updating the 
K-term from the /Lt-terms, fj, < k 7 instead of computing it from scratch. 

We first express ™Fq a ^ (oi :p ; bi :q ; X) in terms of the Jack function J« (X), which 
is normalized so that the coefficient of X1X2 ■ ■ ■ x\ K \ in Jk q) (X) — jfc*\xi,X2, ■ ■ ■ , x n ) 
equals {\k\)\ [23 Thm. 1.1]. The Jack functions C K a) (X) and J K a) (X) are related 
as: 

(3.1) gWffl= ~ ff ^'W, 
where 

(3.2) Jk = n K(i,j)h* K (i,j), 

and h*.(i,j) = k'j — i + a{ni — j + 1) and h%(i,j) = Kj — i + 1 + a(/«j — j) are the 
upper and lower hook lengths at (i,j) € n, respectively. 
Denote 
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Since J^ a \X) = when K n +i > 0, we need to sum only over partitions k with at 
most n parts: 

(3.4) ™F^(a 1:p ;b 1:q] X) = ]T Q,J ( K a) {X). 

|^|<m, K n ^\— 

When computing 1(3.4(1 . we recursively generate all partitions |k| < m in such a 
way that consecutively generated partitions differ in only one part. Therefore it is 
convenient to introduce the notation 

— (^lj ■ ■ • i K i— 1? K i 1? ■ • ■) 

for any partition k such that Ki > /tj+i. 

In the following subsections we derive formulas for updating the /c-term in (13.4(1 
from the /z-terms, fi < n. 

3.1. Updating the Coefficients Q K . We update Q K from Q K{i) using the follow- 
ing lemma. 

Lemma 3.1. 

,„ Q K _ IljLiCgj +c) rr 1 (flj - a)ej ^ Zj - /j 

1 ' J Q« w nf=i + c) /J 5, (e, + a) 1 = 1 I, + ft, ' 

where c = + Ki — 1, d = — i, ej = d — ja + Kj, gj = ej + 1, /j = 

/•Cj-a — j — d, hj = fj + a, and ij = /ij/j . 

Proof. From JOJ, (a)L Q) = (a)&^ ■ (a - (i - l)/a + - 1), and from EHll . 

which along with ((3.3(1 imply 1(3.5(1 . □ 

The seemingly complicated notation of Lemma l3.1l is needed in order to minimize 
the number of arithmetic operations needed to update Q K . A straightforward evalu- 
ation of 1(3.3(1 costs 6|«|(2+p+g); in contrast, 1(3.5(1 costs only 2(p+q) + lOKi + 9i — 11 
arithmetic operations. 

3.2. Updating the Jack Function. When k = (0), J^(xi, . . . ,x n ) — 1. For 

K > (0), we update Jk(x%, . . . ,x n ) from j\T\xi, . . . ,x r ), fi < k, r < n. 

When X is a multiple of the identity we have an easy special case [201 Thm. 5.4]: 

J^(xl)=x^ J] („-(i-i)+ a (j-l)). 

Therefore we can update (xl) from J/fcl (a;/) as 

(3.7) J^(xl) = J<«> (xi) ■ x • (n - i + 1 + a( Ki - 1)). 

In the general case, we update the Jack function using the identity (see, e.g., 
[201 Prop. 4.2]): 

(3.8) J^{ Xl ,x 2 ,...,x n ) = J2jf\x 1 ,X2,...,x n -x)xW^p KI „ 
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where the summation is over all [i < k such that k/ \x is a horizontal strip, and 

( J ^"n^p^j)' Wh6re otherwise. 

The skew partition is a horizontal strip when K\ > /ii > K2 > 7^2 > • ■ ■ j—l 
p. 339]. 

We borrow the idea for updating the Jack function from [3] , but make two impor- 
tant improvements. We only update the coefficients /3 KAt , and store the precomputed 
Jack functions much more efficiently than in PJ. 

The coefficients (3 KfJ , are readily computable using 1)3. 9J1 at the cost of 6(|«| + |//|) 
arithmetic operations. The following lemma allows us to start with @ KK = 1 and 
update /3 Kfi , k -. from /3 KjU at the cost of only 12k + 6/ik — 7 arithmetic operations. 

Lemma 3.2. Let K,fi, and v = /i( fe ) be partitions such that n/fi and k/v are 
horizontal strips, and n' r = jj! r for < r < k — 1 . Then 

(3.10) §* =a A_^.TTtV + a > TT^fV^ 

^ r— 1 r— 1 r— 1 

where a' = a — 1, i = fc — a/ifc, g = t + 1, « r = g — r + an ri v r = t — r + a/i r , and 
w r = fJ,' r —t— ar. 

Proof. Denote Z = 71^. We have v — fi, except for 1//. = — 1. Since and k/i^ 
are horizontal strips, kJ = fi[ and ^ v[ = fjf, — 1. Then 



(3.11) = II 



We transform the first term of l|3.11|l by observing that B* (i, j) = B£ u (i,j), except 
for j = I; 

n^L = TT B nu(r,l) = -rj fegCM) = -pr k - r + 1 + a( Kr - Z) 
RK 11 Rk f r n 11 Wr /I 11 



KjJ. 



?Ji B^(r, I) 1 = 1 Zi*(r, 11 fc - r + a( Kr -l + l) 



To simplify the second term of i|3.11|) . we observe that BgJi,j) = B% v (i,j), except 
for % = fc and j = I. Also n' r = fi' r = v' r for < r < fc— 1, and n' k = fi' k ^ fi' k — 1 = 



Therefore 



W K( k ' r ^ TT ^( r >0 

^-4 /4 - fc + a(/jfe — r + 1) ^-r 1 fc — r + a(/j r — Z + 1) 
11 n' r -k + a(nk - r) ^1 fc - r + a(/j r - Z) 

which yields lj3~TU|) . □ 
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4. Algorithms for Efficient Evaluation of ^F^ a) (ai :p ;bi :q ;X) 

Algorithm 14. II computes p n F q a \ai :p \ b 1 . q ; X) in the easy special case when X is 
a multiple of the identity. Algorithm 14. 21 handles the general case. 

Both algorithms recursively generate all partitions \k\ < m with at most n parts 
by allowing K\ to take all values 1, 2, . . . , m and, independently, Ki,i = 1, 2, . . . , n, 
to take all values 1,2, ... , /tf-j, subject to the restriction \k\ < m. The coefficients 
Q K are updated using 1)3. 5)) . The Jack functions are updated using 1)3. 7)1 or 1)3.8)1 as 
appropriate. 

4.1. The Case when X is a Multiple of the Identity. 

Algorithm 4.1 (Hypergeometric Function, X = xl). The following algorithm 
computes p n F q (ai- p ] bx- q ; xl). The variables x, n, a, s, a, 6, and k are global, 
function s = hgi(m, a, a, 6, n, x) 
s = 1 

summation(l, 1, m) 
function summation(i, z, j) 

for m = 1 : min(Kj_i, j) (defaults to j for i = 1) 

z = z- x-(n — i + 1 + o;(/t, — 1)) ■ T (where T is the right hand side of 1)3. 5(1 ) 
s = s + z 

if (j > Ki) and (i < n) then 

summation(« + 1, z, j — Ki) 
endif 
endf or 

Ki — 

In Algorithm 14.11 the variable z equals the K-term in 1)3.4)1 ; it is updated us- 
ing 1)3.5)1 and 1)3.7)1 . The parameter j in summation equals m — \k\. 

The hypergeometric function of two matrix arguments 1)1.5)1 is in this case 

p F^(a 1:p ;b 1:q] xI;yI) = p ^ Q) (a 1:p ; b 1:q ; xyl). 

4.2. The General Case. We use the identity 1)3.8)1 to update the Jack function, 
but in order to use it efficiently, we need to store and reuse the Jack functions 
computed earlier. 

Therefore we index all partitions \k\ < m with at most n parts (/c n +i = 0) by 
linearizing the m-tree that they form (each node k = (k\, K2, ■ ■ ■ , Kk) has at most 
m children «2, . . . , Kk, Kfe+i), Kk+i = 1, 2, . . . , min(Kfe, m — \k\), k < n). In 
other words, if P mn = |k| < m, K n+ i = 0}, we assign a distinct integer 

index N K 6 {0,1,..., P mn } to every such partition k. We start by assigning the 
indexes 0, 1, . . . , m to partitions with one part: N/q = i, i = 0, 1, 2, . . . , m. Then, 
recursively, once an index N K has been assigned to a partition k = (k%, . . . ,Kk) 
with k parts, we assign m — |k| consecutive unassigned indexes to the partitions 
(ki, . . . , Kk, k/c+i), i^k+i =1)2,..., m— \k\. We record the tree structure in an array 
D such that D(N^ Kl: ,.. yKk )) — Nt Klt .. Kk x). Now given k, we can compute N K by 
starting with X( Kl ) = k\ and using the recurrence 

(4-1) N (Ku ..., Ki) = D(N {Kl _ Kz _ l} ) +Ki-1. 

We store every computed Jk (xi-a) (\k\ < m, K n+ \ = 0, % = 1, 2, . . . ,n) in the 
(N K , i)th entry of an P mn x n array, which we call "J" in Algorithm 14. 21 below. 
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We compute the value of P m n as follows. Let Pk(i) be the number of partitions 
of i with exactly k parts. The Pk(i), i — 1,2, ... , to, k = 1,2,..., min(m, n), are 
computed using the recurrence Pk(i) = Pk-i(i ~ 1) +Pfc(* — k) [211 p. 28]. Then 

m min(m,n) 

(4.2) Pnm=J2 £ 

i=l fe=l 

Next, we present our main algorithm. 

Algorithm 4.2 (Hypergeometric Function). The following algorithm computes 
™Fq (ai ;p ',bi:q',X), where X = diag(a;). The variables m,n,a,b,a,x,J,D,H,K, 
N K , fi, Np, and s are global, 
function s = hg(m, a, a, 6, x) 
Compute P mn using (|4.2|l 

ri = length(.T); H = m + I; s = 1; I? = zeros(P„ m , 1) 

J = zeros(P mn , n); J(l, :) = 1 {Jk*\x x-j) is stored in J(N K ,i)) 

summation(l, 1, 1, to) 

function summation(i, z, j) 

r = min(/Cj_i, j) (defaults to j for i = 1) 

for Kj = 1 : r 

if (Ki = l) and (i > 1) then D(N K ) = H; N K = H; H = H + r 

else 7V K = iV K + 1 

endif 

z = z • T (where T is the right hand side of <|3.5|) ) 

if = 1 then 

J(JV K , l) = xi(l + a(/si - 1)) • J(iV K - 1, 1) 
endif 

for t = 2:n do jack(0, 1, 0, t) (computes (xi-.t)) 

endf or 

s = s + z ■ J(N K , n, 1) 

if (j > Ki) and (i < n) then 

summation(i + 1, z,j — Ki) 
endif 
endf or 

function jack(fc, /3 Kfl , c, t) 
for i = k : fi[ 

if (fc > 0) and (/j; > fn+l) 
d = N^ 

Hi = /j,i — 1; Compute iV M using l|4.1|l 

Update /3 KjU using l|3.10jl 

if fii > then jack(i, /3 re/J , c + 1, t) 

else J(JV„,t) = J(JV B ,t) +/3 KM • J(JV M , i - 1) • x^ 1 

endif 

= ^ + 1; N^ = d 
endif 
endf or 

if fc = then J(7V K , i) = J(N K , t) + J(N K , t - 1) 

else J(N K , t) = J(N K , t) + • 

endif 
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Algorithm 14. 21 generates all partitions \n\ < m with at most n parts analogously 
to Algorithm 14.11 The parameter z in summation is now just Q K . 

For every k the function jack recursively generates all partitions fi < k such 
that n/fi is a horizontal strip. The Jack function is computed using l|3.8[l : /3 KjU is 
updated using (|3.10() . The parameter c equals |k| — at all times. 

4.3. Implementation notes. In our implementation of Algorithms 14.11 and 14.21 

in |13) we made a few fairly straightforward, but important improvements: 

(1) we precompute the values of xj for i = 1, 2, . . . , n, and j = 1, 2, . . . , m; 

(2) we keep and update the conjugate partitions k' and yl along with k and /x, 
so we never need to recover k 1 and // when computing (|3.5|) and (|3.1U|) : 

(3) when two sets of arguments (x\. n and yi :n ) are passed, the hypergeometric 
function of two matrix arguments, (|1.5fl . is computed; 

(4) if a vector x — (ti,...,t r ) is passed as a parameter in Algorithm 14.11 
the output is also a vector with the values of ™Fq a \ai; P ;bi :q ;tjI) for j — 
1,2,. ...r. 

5. Complexity Analysis 

Algorithm 14.11 (the case X = xl) costs 0(P mn ) arithmetic operations (where 
again, P mn is the number of partitions k, |k| < m with at most n parts). 

To bound the complexity of Algorithm 14.21 (general case) we observe that the 
formula 1(3. 8|) represents the summation of at most P mn terms (in fact a lot less, 
but we have been unable to obtain a better bound than P mn that is easy to work 
with). We use (|3.8|l to compute the Jack function J« (xi-.t) for all |k| < m and 
t < n, i.e., P m n • n Jack functions, each of which costs at most 0(P mn ) arithmetic 
operations. The overall cost of Algorithm 14. 21 is thus bounded by 

0(P mn ■ n). 

There is no explicit formula for P mn , so we use Ramanujan's asymptotic for- 
mula (111 p. 116] for number of partitions of m 

P(m) - ^^/g exp (nV 2m /3) 

to obtain 

rn 

P mn < J2p( 1 ) ~ ( ex P (*V2m/3)) , and P mn - O (exp (2Tr^/2m/3)^ . 

i=l 

Therefore the complexity of Algorithm 14. 21 is linear in n and subexponential in m. 

6. Numerical Experiments 

We performed extensive numerical tests to verify the correctness and complex- 
ity of our Algorithms 14.11 and 14.21 We compared the output of our algorithms for 
qFq and i-Fj with explicit expressions ( subsection 16.1(1 . We also compared the 
probability distributions of the eigenvalues of certain random matrices (which are 
expressed as hypergeometric functions) against the results of Monte-Carlo experi- 
ments fsubsection l6.2|l . Finally, we present performance results in subsection 16.31 
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6.1. Explicit Expressions. We compared the output of Algorithms 14.11 and 14.21 

against the expressions p. 262]: 

(6.1) oFt\X) = e^ x ; 

(6.2) iFk a) (a\X) = det(I-X)- ( \ 

for n = 10 and random uniformly distributed values Xi G [0, 4], i = 1, 2, . . . , n. For 
m = 30 it took less than one second per test. The results agreed to at least 12 
decimal digits with Qfi.ljl . and at least 8 decimal digits with ()6.2J) . reflecting the 
slower convergence of 1)6.2(1 . 

6.2. Eigenvalue Statistics. We tested Algorithms 14. II and 14.21 against the eigen- 
value statistics of the /3-Laguerre and Wishart matrices. 

The n x n [3-Laguerre matrix of parameter a is defined as L = BB T , where 



a > !)■ 

The n x n Wishart matrix with I degrees of freedom (I > n) and covariance 
matrix E, A = W n (l, E), is defined as A = E 1 / 2 ■ Z T ■ Z ■ E 1 / 2 , where = iV(0, 1) 
for i = 1,2, j = 1,2,..., n. 

The eigenvalue distributions of A and L are the same when (3 — 1 or 2, a = |, 
and E = /. 

The cumulative distribution functions of the largest eigenvalues, \l and A ,4, of 
L and A, are 

P^<x) = e r (a VJ + 1) @ x^(a;a + ^ + l;-M. 
^<x) = " i + 2 J det^E-^^^d;^;-^- 1 ), 

ln \ 2 J 

respectively gj Thm. 10.2.1, p. 147], 17, Thm. 9.7.1, p. 420], where a = 2/(3, and 
r„ is the multivariate Gamma function of parameter a: 

(6.3) T^{c)=Tr^f\r(c-—) for5i( c )>— . 

±A - \ a ) a 

We use the Kummer relation \I7\ Thm. 7.4.3, p. 265] 

1 F^(a; c; X) = e tlX ■ ^(c - a; c; -X) 

to obtain the equivalent, but numerically more stable expressions 

^ <x)= ; \ a J/^ (| e"* x > (^1 + 1; a + + 1; M ; 

p(2)/n+l\ 

P(A^ < x) = " l + 2 +1 j det (I,E-^) i/2 eM-f ^) lif ) (2 ±i ; a±m ; i^-i) , 

1 " I 2 ) 

which we plot on Figure ^ along with the Monte-Carlo results from a sample of 
10000 random matrices. 



B = 



Xj3(n-l) X.2a-fi 



Xp X2o-fl(n-l) 
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cdf of X of 3x3 2-Laguerre, a=3 
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Figure 1. The c.d.f. of A max of the /3-Laguerre matrix (top plots) 
and the Wishart matrix (bottom plots). 



Next we consider the distribution of the smallest eigenvalue of L and that of 
txA. 

Ifc = a — %{n — 1) — 1 is a nonnegative integer, then the probability density 
function of the smallest eigenvalue of L is (see, e.g., Thm. 10.1.1, p. 146]): 



(6.4) 



T?(2/0)( „ m 
2-ffJ {-C,P- 



-Jn-l) 



Since — c is a nonpositive integer, the series expansion of a-fn m jEH terminates, 
even though it diverges in general. 

The probability density function of tr A is 

gM +ki2X (u) 



/(«) 



det(A" 1 E)-'/ 2 ^ 



k=0 

oo 



k) 



■sl-fc 



(6.5) = det(A- 1 S)- i / 2 E %+MA (u)E 2 fc J 7 1 ■ J« 2) (/ - AS" 1 ), 
where 



k=0 Khfe 
p -u/2X u r-l 



(u) 



(« > 0), 



(2A)T(r) 

and A is arbitrary H3 P- 341]. We follow the suggestion by Muirhead to use 
A = 2AiA;/(Ai + A/), where Ai and A; are the largest and smallest eigenvalues of 
E, respectively. Although the expression l|6.5JI is not a hypergeometric function of 
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pdf of X . of 5x5 0.5-Laguerre, a=5 
r mm 3 



pdf of X . of 5x5 6-Laguerre, a=1 6 
r mm 3 




pdf of tr(A), A=W 3 (4,diag(0.25,0.5,1)) 



0.15 




0.05 




pdf of tr(A), A=W 4 (5,diag(0.5,1 ,1 .1 ,2)) 




80 100 



Figure 2. The p.d.f. of A m i n of the /3-Laguerre matrix (top plots) 
and the p.d.f. of the trace of the Wishart matrix (bottom plots). 



a matrix argument, its truncation for \k\ < m has the form l|1.6|l . and is computed 
analogously. 

We plot (|6.4|) and (|6.5[) on Figure [3 and compare the theoretical predictions of 
these formulas with experimental data. 



6.3. Performance Results. In Figure EH we demonstrate the efficiency of Algo- 
rithms ^3 an d ^21 on an 1.8GHz Intel Pentium 4 machine. 

In the left plot we present the performance data for Algorithm 14.11 f whose com- 
plexity is independent of the size n of the matrix X — xl) . Its efficiency is evident — 
we need to sum beyond partitions of size m = 52 before this algorithm takes a full 
second (for reference, the fc! in the denominator of the K-term in then reaches 
up to 52! fa 8 • 10 67 ). 

Algorithm 14.21 is also very efficient. The right plot of Figure [3| demonstrates 
clearly its linear complexity in n. It also takes at most a few seconds on matrices 
of size n < 120 and partitions of size m < 30. 



7. Conclusions and Open Problems 



We have presented new algorithms for computing the truncation of the hyperge- 
ometric function of a matrix argument. They exploit the combinatorial properties 
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Time to compute cT^o ) Time to compute ™^o 2) (X) 




20 30 40 50 O 50 100 

max size partition — m n = size of X 



Figure 3. The performance of Algorithms I4.ll (left) and l4.2l (right). 

of the Pochhammer symbol and the Jack function to achieve remarkable efficiency, 
and have lead to new results H>] . 

Several problems remain open, among them automatic detection of convergence. 
The K-term in l|3.4|) does approach zero as \k\ — > oo, but it need not monotonically 
decrease. Although we have 

TT (ffj ~ a ) e j TT l i - fj < 1 

in (|3.5|l . it is not always true that Q k Jk < Qku^Jk^i and it is unclear how to tell 
when convergence sets in. 

Another open problem is to determine the best way to truncate the series <|l.ll) . 
Our choice to truncate it for |k| < m seems to work well in practice, but one can 
imagine selecting a partition A and truncating for k < A instead of, or in addition 
to | AC | < TO. 
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